## "Public Preferences for International Law Compliance: 
##  Respecting Legal Obligations or Conforming to Common Practices?"
##
##
## Reproduce: 
##      - Figure 2 
##
## Input:  
##      - data/data_cleaned.csv
##
## Output: 
##      - f2.pdf
##
## Saki Kuzushima
## 01/06/2023

source("code/functions_f2.R")
data <- read.csv(file="data/data_cleaned.csv")

# output table
item <- c("surname3", 
          "whale1", "hate1", "death1")
res <- expand.grid(item=item, stringsAsFactors=F)
res$lower <- 0
res$point <- 0
res$upper <- 0
res$p <- 0
res$se <- 0

for (i in seq(1, nrow(res))){
  item <- res[i, "item"]
  res[i, c("lower", "point", "upper", "p", "se")] <- 
    diffmean_block(data=data,  item=item)
}


# Multiple testing correction
res$p_bh <- p.adjust(res$p, method="BH")

# adaptive shrinkage
ash_fit <- ash(betahat=res$point, 
               sebetahat=res$se, 
               prior="uniform", 
               mixcompdist="unif",
               grange=c(-5, 5),
               outputlevel=2)
ash_result <- ash_fit$result
ash_ci <- ashci(ash_fit)
res_ash <- data.frame(item=res$item,
                      lower=ash_ci[,1], 
                      point=ash_result[,"PosteriorMean"],
                      upper=ash_ci[,2])
# Figure
x_index <- c(1.25, 1.75, 3, 4.25, 4.75, 6, 7.25, 7.75, 9, 10.25, 10.75) 
point <- c(res$point[1], res_ash$point[1], NA, 
           res$point[2], res_ash$point[2], NA, 
           res$point[3], res_ash$point[3], NA, 
           res$point[4], res_ash$point[4])
upper <- c(res$upper[1], res_ash$upper[1], NA, 
           res$upper[2], res_ash$upper[2], NA, 
           res$upper[3], res_ash$upper[3], NA, 
           res$upper[4], res_ash$upper[4])
lower <- c(res$lower[1], res_ash$lower[1], NA, 
           res$lower[2], res_ash$lower[2], NA, 
           res$lower[3], res_ash$lower[3], NA, 
           res$lower[4], res_ash$lower[4])


# Style memo:
# Point style: hollow for uncorrected, filled for corrected ones
# Note that pch=18 (filled diamond) looks smaller than others. 
# So the cex for pch=18 is set to 2.0; others 1.5
# CI color: grey for non-significant ones; Black for significance
point_style <- c(1,16,NA,2,17,NA,0,15,NA,5,18)
pdf("results/main/f2.pdf")
par(mfrow=c(3,1), mar=c(1,4,1,1)+0.1, oma=c(2,1,3,1)+0.1)
plot(x=x_index, y=point, type="n",
     xlim=c(0.1, 11.9), ylim=c(-0.5, 0.5),
     xaxs="i", yaxs="i",
     cex.lab=1.25,
     xaxt="n", yaxt="n",
     xlab="", ylab="Constitution - Int'l law ",
     pch=point_style,
     cex=1.5,
     col="gray50")
abline(h=0, col="red")
abline(v=3, lty=3)
abline(v=6, lty=3)
abline(v=9, lty=3)
axis(side=3, at=c(1.5, 4.5, 7.5, 10.5), labels=c("Surname", "Whale", "Hate Speech", "Death Penalty"),
     tick=F,  line=0,  outer=T, cex.axis=1.25) 
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(1,4,7,10)], 
       y=point[c(1,4,7,10)],
       pch=point_style[c(1,4,7,10)],
       cex=1.5,
       col="grey50") 
points(x=x_index[c(2,5,8,11)], 
       y=point[c(2,5,8,11)],
       pch=point_style[c(2,5,8,11)],
       cex=c(1.5, 1.5, 1.5, 2),
       col="grey50") 
points(x=x_index[10], 
       y=point[10],
       pch=point_style[10],
       cex=1.5,
       col="black") 
segments(x0=x_index[c(1,4,7,10)], 
         x1=x_index[c(1,4,7,10)],
         y0=lower[c(1,4,7,10)],
         y1=upper[c(1,4,7,10)], col="grey50") 
segments(x0=x_index[c(2,5,8,11)], 
         x1=x_index[c(2,5,8,11)],
         y0=lower[c(2,5,8,11)],
         y1=upper[c(2,5,8,11)], col="grey50") 
segments(x0=x_index[10], 
         x1=x_index[10],
         y0=lower[10],
         y1=upper[10], col="black") 
dev.off()
